#include "../include/Equation.h"
#include "../include/Interp.h"
#include <iostream>

using namespace std;

void splitline()
{
    cout << "---------------------------------------" << endl;
}

double f1(double x)
{
    return 1.0 / (1 + x * x);
}

double f2(double x)
{
    return 1.0 / (1 + 25 * x * x);
}

void exerciseA()
{
    cout << "Test output" << endl;
    Polynomial poly1({1, 2, 3, 4, 5, 6, 7, 8, 9, 10});
    Polynomial poly2({1, 2, 3, 4, 5});
    cout << "poly1 = " << poly1 << endl;
    cout << "poly2 = " << poly2 << endl;

    cout << "Test addition, subtraction, multiplication, evaluation" << endl;
    cout << "poly1 + poly2 = " << poly1 + poly2 << endl;
    cout << "poly2 - poly1 = " << poly2 - poly1 << endl;
    cout << "poly1 * poly2 = " << poly1 * poly2 << endl;
    cout << "poly1(1) = " << poly1.caculate(1) << endl;
    cout << "poly1(2) = " << poly1.caculate(2) << endl;

    cout << endl
         << "Test derivation and integration" << endl;
    cout << "p(x) = " << poly1 << endl;
    cout << "p'(x) = " << poly1.derivation() << endl;
    cout << "P(x) = " << poly1.integration() << endl;

    cout << endl
         << "Test Hermite Interpolation" << endl;
    INTERPC interpC;

    interpC.points = {0, 1, 1, 1};
    interpC.interpSites[0] = {0};
    interpC.interpSites[1] = {1, 7, 42};
    NewtonInterp newtonInterp;

    cout << "Overwrite: " << endl;
    newtonInterp.overwriteInterp(interpC);
    cout << newtonInterp << endl;

    interpC.points = {2, 2};
    interpC.interpSites[2] = {128, 448};
    newtonInterp.incrementInterp(interpC);
    cout << "Incremently: " << endl;
    cout << newtonInterp << endl;

    cout << "Test Neville_Aitken Algorithm" << endl;
    interpC.points = {0, 1, 2, 3, 4};
    interpC.interpSites[0] = {6};
    interpC.interpSites[1] = {-3};
    interpC.interpSites[2] = {-6};
    interpC.interpSites[3] = {9};
    interpC.interpSites[4] = {15};
    cout << NewtonInterp::Neville_Aitken(interpC) << endl;
}

void exerciseB(int n)
{
    double x = -5;
    INTERPC interpC;
    for (int i = 0; i < n + 1; i++)
    {
        interpC.points.push_back(x + 10.0 * i / n);
        interpC.interpSites[interpC.points[i]] = {f1(interpC.points[i])};
    }
    NewtonInterp newtonInterp;
    newtonInterp.overwriteInterp(interpC);
    cout << newtonInterp << endl;
}

void exerciseC(int n)
{
    NewtonInterp newtonInterp;
    newtonInterp.Chebyshev(n, f2);
    cout << newtonInterp << endl;
}

void exerciseD()
{
    vector<double> time = {0, 3, 5, 8, 13};
    vector<double> distance = {0, 225, 383, 623, 993};
    vector<double> speed = {75, 77, 80, 74, 72};

    INTERPC interpC1;
    interpC1.points = time;
    for (int i = 0; i < (int)time.size(); i++)
    {
        interpC1.interpSites[time[i]] = {distance[i]};
    }

    cout << "Distance:" << endl;
    NewtonInterp newtonInterp1;
    newtonInterp1.overwriteInterp(interpC1);
    cout << newtonInterp1 << endl;
    cout << "Distance for t = 10s: " << newtonInterp1.get_poly().caculate(10) << endl;

    INTERPC interpC2;
    interpC2.points = time;
    for (int i = 0; i < (int)time.size(); i++)
    {
        interpC2.interpSites[time[i]] = {speed[i]};
    }

    cout << "Speed:" << endl;
    NewtonInterp newtonInterp2;
    newtonInterp2.overwriteInterp(interpC2);
    cout << newtonInterp2 << endl;
    cout << "Speed for t = 10s: " << newtonInterp2.get_poly().caculate(10) << endl;

    Polynomial sp = newtonInterp2.get_poly().derivation();
    cout << endl
         << "Derivation of speed:" << endl;
    cout << sp << endl;

    //speed(x) = -3.2891 + 4.47718x^1 + -1.05731x^2 + 0.0610256x^3
    double (*f)(double) = [](double x) -> double
    {
        return -3.2891 + 4.47718 * x - 1.05731 * x * x + 0.0610256 * x * x * x;
    };

    Bisection bisection(f);
    double point = bisection.solve(3, 8);
    cout << "Zero position: " << point << endl;
    cout << "Max value: " << newtonInterp2.get_poly().caculate(point) << endl;
}

void exerciseE()
{
    vector<double> day = {0, 6, 10, 13, 17, 20, 28};
    vector<double> sp1 = {6.67, 17.3, 42.7, 37.3, 30.1, 29.3, 28.7};
    vector<double> sp2 = {6.67, 16.1, 18.9, 15.0, 10.6, 9.44, 8.89};

    INTERPC interpC1;
    interpC1.points = day;
    for (int i = 0; i < (int)day.size(); i++)
    {
        interpC1.interpSites[day[i]] = {sp1[i]};
    }

    cout << "Sp1:" << endl;
    NewtonInterp newtonInterp1;
    newtonInterp1.overwriteInterp(interpC1);
    cout << newtonInterp1 << endl;

    INTERPC interpC2;
    interpC2.points = day;
    for (int i = 0; i < (int)day.size(); i++)
    {
        interpC2.interpSites[day[i]] = {sp2[i]};
    }

    cout << "Sp2:" << endl;
    NewtonInterp newtonInterp2;
    newtonInterp2.overwriteInterp(interpC2);
    cout << newtonInterp2 << endl;

    cout << "After 15 Days:" << endl;
    cout << "Sp1: " << newtonInterp1.get_poly().caculate(43) << endl;
    cout << "Sp2: " << newtonInterp2.get_poly().caculate(43) << endl;
}

int main(int argc, char *argv[])
{
    cout << "ExerciseA" << endl;
    exerciseA();
    splitline();

    cout << "ExcercseB" << endl;
    for (int i = 0; i < 4; i++)
    {
        cout << "n = " << 2 * (i + 1) << endl;
        exerciseB(2 * (i + 1));
    }
    splitline();

    cout << "ExerciseC" << endl;
    for (int i = 0; i < 4; i++)
    {
        cout << "n = " << 5 * (i + 1) << endl;
        exerciseC(5 * (i + 1));
    }
    splitline();

    cout << "ExerciseD" << endl;
    exerciseD();
    splitline();

    cout << "ExerciseE" << endl;
    exerciseE();

    return 0;
}
